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Abstract 


Lagrangian  Stochastic  (LS)  particle  models  have  proven  to  be  a  useful  computational  tool  for  the 
description  and  prediction  of  dispersion  of  pollutant  releases  in  complex  meteorological  situations 
(e.g.,  space-  and  time-varying  situations  pertaining  to  complex  flow  and  turbulence).  However, 
simulating  the  emitted  pollutant  by  following  the  trajectories  of  many  “marked”  fluid  elements 
released  from  the  source  distribution  brings  up  the  difficulty  of  the  correct  estimation  of  the  mean 
concentration  of  the  dispersing  pollutant  from  the  particle  trajectory  information.  Recently,  the 
density  kernel  estimation  method  has  been  proposed  and  applied  successfully  to  estimate  mean 
concentrations  from  Lagrangian  Stochastic  particle  models.  However,  the  computational  effort 
needed  by  this  method  increases  as  TV2  (assuming  the  number  of  receptor  locations  Nr  at  which  the 
concentration  is  required  is  comparable  to  the  number  of  fluid  particles  Np  used  in  the  trajectory 
simulation,  soNr~Np~N)  and,  in  consequence,  the  method  has  not  been  widely  used  because  of 
the  significant  computer  resources  required.  Here,  we  describe  a  novel  algorithm  for  calculating 
the  kernel  estimate  of  the  mean  concentration  field  whose  computational  complexity  scales  only  as 
N.  The  technique  uses  a  tesselation  (subdivision)  of  space  in  cubic  cells  of  side  length  h  (where  h 
is  the  bandwidth  of  the  kernel  function),  and  then  associates  a  linked-list  data  structure  with  each 
cell  that  is  used  as  a  bookkeeping  device  to  keep  track  of  the  “marked”  fluid  particles  in  that  cell. 
The  fast  approach  developed  here  has  been  verified  by  comparing  results  with  the  direct 
implementation  of  the  kernel  estimator  and  with  the  conventional  box-counting  estimator  for  the 
mean  concentration  field. 


Resume 


On  a  prouve  que  les  modeles  de  particules  langrangiens  stochastiques  (LS)  sont  des  outils  utiles  de 
calcul  pour  decrire  et  predire  la  dispersion  des  emissions  de  polluants  dans  des  situations 
meteorologiques  complexes  (par  ex.  :  des  situations  dans  l’espace  et  le  temps  qui  varient  selon  des 
flux  et  des  turbulences  complexes).  Simuler  le  polluant  emis  en  suivant  les  trajectoires  de 
beaucoup  d’elements  de  fluides  «  marques  »  qui  ont  ete  emis  a  la  repartition  des  sources, 
augmente  cependant  la  difficulty  d’estimer  correctement  la  concentration  moyenne  du  polluant 
disperse,  a  partir  de  I’information  sur  la  trajectoire  de  la  particule.  La  methode  d’estimation  a 
partir  du  noyau  de  densite  a  ete  recemment  proposee  et  appliquee  avec  succes  pour  estimer  les 
concentrations  moyennes  des  modeles  de  particules  langrangiens  stochastiques.  Cependant, 

1’ effort  de  calcul  requis  par  cette  methode  augmente  par  un  facteur  N2  (dans  l’hypothese  oil  le 
nombre  de  lieux  recepteurs  Nr  auxquels  la  concentration  est  requise  est  comparable  au  nombre  de 
particules  de  fluides  Nr  t utilise  dans  la  simulation  de  la  trajectoire,  tel  que  Nr »  Np  ~  N)  et,  par 
consequent,  cette  methode,  exigeant  des  ressources  informatiques  tres  importantes,  n’a  pas  ete  tres 
utilisee.  Nous  decrivons  ici,  un  nouvel  algorithme  qui  calcule  l’estimation  du  noyau  de.densite  du 
champ  de  la  concentration  moyenne  et  dont  la  complexity  des  calculs  est  reduite  a  N.  La  technique 
utilise  une  tessellation  (subdivision)  de  1’espace  en  cellules  cubiques  dont  la  longueur  des  cotes  h 
(h  etant  la  largeur  de  bande  de  la  fonction  du  noyau)  et  puis  associe  une  structure  de  donnees  en 
chaine  avec  chaque  cellule;  cette  structure,  etant  utilisee  comme  engin  de  comptabilite,  suit  les 
particules  de  fluides  «  marquees  »  dans  cette  cellule.  Cette  methode  rapide,  mise  au  point  ici,  a  ete 
verifiee  en  comparant  des  resultats  avec  1’ implementation  directe  de  l’estimateur  de  noyau  et  avec 
1‘estimateur  classique  de  «  comptage  de  boites  »  pour  le  champ  de  concentration  moyenne. 
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Executive  summary 

Background 

There  are,  quite  rightly,  growing  concerns  world-wide  about  the  dangers,  both  actual  and 
potential  of  the  use  (release)  of  chemical  and  biological  warfare  (CBW)  agents  into  the 
atmosphere.  The  increased  awareness  and  importance  accorded  by  the  public  world-wide  and 
their  governments  to  maintain  appropriate  defences  against  CBW  agents  in  the  environment, 
the  prediction  of  casualties  and  human  performance  degradation  resulting  from  such  releases, 
and  the  development  of  operational  procedures  and  regulations  to  control,  mitigate,  and 
monitor  the  fate  of  CBW  agents  in  the  atmosphere  will  require  mathematical  modelling  of 
atmospheric  dispersion  of  these  agents.  This  modelling  is  necessary  in  order  to  answer 
political  and  planning  questions  such  as:  what  is  the  extent  of  the  hazard  zone  corresponding 
to  the  release  of  a  CBW  agent  and  what  is  the  effect  of  this  release  on  exposed  personnel 
within  the  hazard  region;  and  what  contingency  measures  would  be  needed  to  deal  with  a 
potential  release  of  a  CBW  agent.  Accurate  prediction  of  the  dispersion  of  contaminants 
released  into  the  complex  atmospheric  boundary  layer  is  a  particularly  challenging  problem. 
Compared  to  conventional  techniques  for  the  calculation  of  dispersion  (e.g.,  similarity  theory, 
statistical  theory,  or  eddy-diffusivity  methods),  Lagrangian  Stochastic  (LS)  models  of  particle 
motions  have  proven  to  be  a  successful  and  flexible  tool  in  the  description  and  prediction  of 
contaminant  dispersion  in  complex  meteorological  situations  (incorporating  in  a  natural 
manner  the  effects  of  inhomogeneities,  unsteadiness  and/or  non-Gaussianity  in  turbulent 
velocity  distributions  associated  with  the  necessarily  complex  flow  and  turbulence 
encountered  in  the  real-world  atmosphere). 

The  major  disadvantage  of  the  LS  modelling  approach  is  the  need  for  large  computational 
resources  to  run  the  model  for  typical  applications  in  order  to  obtain  statistically  meaningful 
solutions  of  contaminant  dispersion  and  mean  concentrations  for  sources  in  typical  CBW 
agent  releases.  To  facilitate  the  application  of  state-of-the-art  LS  models  for  CBW  agent 
dispersion  predictions,  we  developed  a  computationally  efficient  algorithm  for  mean 
concentration  estimation  based  on  the  density  kernel  method  that  can  be  used  in  conjunction 
with  LS  models  of  dispersion. 

Principal  results 

Lagrangian  Stochastic  particle  models  have  proven  to  be  a  useful  computational  tool  for  the 
description  and  prediction  of  dispersion  of  pollutant  releases  in  complex  meteorological 
situations  (e.g.,  space-  and  time-varying  situations  pertaining  to  complex  flow  and 
turbulence).  However,  simulating  the  emitted  pollutant  by  following  the  trajectories  of  many 
“marked”  fluid  elements  released  from  the  source  distribution  brings  up  the  difficulty  of  the 
correct  estimation  of  the  mean  concentration  of  the  dispersing  pollutant  from  the  particle 
trajectory  information.  Recently,  the  density  kernel  estimation  method  has  been  proposed  and 
applied  successfully  to  estimate  mean  concentrations  from  LS  particle  models.  However,  the 
computational  effort  needed  by  this  method  increases  as  N2  (assuming  the  number  of  receptor 
locations  Nr  at  which  the  concentration  is  required  is  comparable  to  the  number  of  fluid 
particles  Np  used  in  the  trajectory  simulation,  so  Nr~Np~N)  and,  in  consequence,  the  method 
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has  not  been  widely  used  because  of  the  significant  computer  resources  required.  Here,  we 
describe  a  novel  algorithm  for  calculating  the  kernel  estimate  of  the  mean  concentration  field 
whose  computational  complexity  scales  only  as  N.  The  technique  uses  a  tesselation 
(subdivision)  of  space  in  cubic  cells  of  side  length  h  (where  h  is  the  bandwidth  of  the  kernel 
function),  and  then  associates  a  linked-list  data  structure  with  each  cell  that  is  used  as  a 
bookkeeping  device  to  keep  track  of  the  "marked”  fluid  particles  in  that  cell.  The  fast 
approach  developed  here  has  been  verified  by  comparing  results  with  the  direct 
implementation  of  the  kernel  estimator  and  with  the  conventional  box-counting  estimator  for 
the  mean  concentration  field. 

Significance  of  results 

The  development  and  implementation  of  a  fast  linked  list  based  algorithm  for  a  density  kernel 
estimate  of  the  mean  concentration  field  has  been  described,  allowing  the  user  to  calculate  a 
set  of  kernel  estimates  of  concentration  at  Nr  receptor  locations  in  order  0(Nr)  work  as 

opposed  to  order  0(N2r )  effort  (assuming  the  typical  situation  where  the  number  of  fluid 
particles  Np  used  in  the  trajectory-simulation  satisfies  Np  >Nr).  The  algorithm  makes  feasible 
the  use  of  LS  models  on  a  personal  computer  for  the  routine  calculation  of  dispersion 
associated  with  complex  CBW  agent  release  scenarios  into  the  atmosphere. 

Future  work 

There  are  several  ways  in  which  the  code  for  the  fast  linked  list  based  algorithm  for  the 
density  kernel  estimator  of  the  mean  concentration  field  can  be  generalized  and/or  made  more 
efficient.  The  linked  list  data  structure  used  here  is  appropriate  for  the  case  where  the 
bandwidth  of  the  density  kernel  associated  with  each  particle  is  fixed.  In  principle,  the 
bandwidth  h  need  not  be  constant  and  for  many  applications  it  may  be  useful  to  allow  a 
variable  smoothing  length  that  is  dynamically  adapted  so  that  the  number  of  neighboring  fluid 
particles  within  the  support  of  a  kernel  function  centered  on  any  particular  particle  remains 
constant.  The  generalization  of  the  proposed  fast  kernel  estimator  for  variable  kernel 
bandwidth  h  will  probably  require  that  the  linked  list  data  structure  used  here  be  replaced  by  a 
hierarchy  tree  data  structure  that  can  be  adopted  to  suit  the  needs  of  a  variable  bandwidth.  In 
the  application  of  the  fast  kernel  estimator,  the  code  can  be  made  even  more  efficient  by  pre- 
computing  the  values  of  the  kernel  function  for  a  large  number  of  representative  points  and 
storing  the  results  in  a  pre-computed  kernel  function  value  look-up  table.  Finally,  it  may  well 
be  worth  exploring  the  parallelization  of  the  code  so  that  it  can  be  executed  on  computers  with 
highly  parallel  architectures  (e.g.,  Beowulf  clusters). 


Shao,  Y.  and  Yee,  E.  (2004).  The  Rapid  Evaluation  of  Mean  Concentration  Fields  in  Lagrangian 
Stochastic  Modelling  Using  a  Density  Kernel  Estimator,  (DRDC  Suffield  TR  2004-1 86), 
Defence  R&D  Canada  -  Suffield. 
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Sommaire 


Contexte 

II  existe,  a  juste  titre,  des  inquietudes  croissantes  au  niveau  mondial,  en  ce  qui  concerne  Ies 
dangers,  a  la  fois  actuels  et  potentiels,  d’utilisation  (emission)  d’agents  chimiques  et  biologiques 
de  guerre  (CBW),  dans  Patmosphere.  Cette  conscience  accrue  et  l’importance  accordee  par  le 
public  et  par  les  gouvemements,  au  niveau  mondial,  le  maintien  des  defenses  appropriees  contre 
les  agents  CBW  dans  l’environnement,  la  prevision  de  victimes,  la  degradation  des  performances 
humaines  resultant  de  telles  emissions  et  l’elaboration  de  procedures  operationnelles  et  de 
reglements  pour  controler,  limiter  et  surveiller  le  sort  des  agents  CBW  dans  1’ atmosphere, 
exigeront  la  moderation  mathematique  de  la  dispersion  de  ces  agents  dans  Patmosphere.  Cette 
modelisation  est  necessaire  pour  faire  face  aux  problemes  politiques  et  de  planification  tels  que  : 
l’etendue  de  la  zone  dangereuse  resultant  de  remission  d’un  agent  CBW  et  l’effet  d’une  telle 
emission  sur  le  personnel  qui  a  ete  expose  a  Pinterieur  de  la  region  dangereuse  ;  et  les  mesures 
d’urgence  requises  pour  gerer  remission  potentieile  d’un  agent  CBW.  Predire  avec  exactitude  la 
dispersion  de  ces  contaminants  emis  dans  le  couche  limite  atmospherique  est  tres  problematique. 
Compares  aux  techniques  classiques  de  calcul  de  la  dispersion  (par  ex. :  theorie  de  la  similitude, 
theorie  statistique  ou  les  methodes  de  diffusivite  des  remous),  les  mouvements  des  modeles  de 
particules  lagrangiens  stochastiques  (LS)  se  sont  prouves  etre  des  outils  souples  qui  reussissent  a 
decrire  et  a  predire  la  dispersion  de  contaminants  dans  des  situations  meteorologiques  complexes 
(incorporant  de  maniere  naturelle  les  effets  d’elements  non  homogenes,  le  manque  de  regularity  et 
/ou  les  distributions  non  gaussiennes  de  vitesse  limite  de  regime  turbulent  associees  avec  les  flux 
et  les  turbulences  complexes  qui  sont  necessairement  rencontrees  dans  la  realite  de  Patmosphere). 

L’ inconvenient  principal  de  cette  methode  de  moderation  LS  est  que  cette  derniere  necessite  des 
ressources  tres  importantes  en  puissance  de  calcul  pour  que  ce  modele  produise  des  solutions  aux 
applications  caracteristiques  qui  soient  statistiquement  valables  en  ce  qui  concerne  la  dispersion  de 
contaminants  et  les  concentrations  moyennes  des  sources,  lors  d’emissions  caracteristiques 
d’agents  CBW.  Un  algorithme  de  traitement  efficace  estimant  la  concentration  moyenne  basee  sur 
la  methode  de  la  densite  du  noyau  a  ete  developpe  pour  faciliter  l’application  des  ces  modeles  LS, 
les  plus  recents  de  la  technique,  aux  predictions  de  dispersion  d’agents  CBW  et  peuvent  etre 
utilises  en  conjonction  avec  les  modeles  LS  de  dispersion. 

Les  resultats  principaux 

On  a  prouve  que  les  modeles  de  particules  langrangiens  stochastiques  sont  des  outils  utiles  de 
calcul  pour  decrire  et  predire  la  dispersion  des  emissions  de  polluants  dans  des  situations 
meteorologiques  complexes  (par  ex.  :  des  situations  dans  Pespace  et  le  temps  qui  varient  selon  des 
flux  et  des  turbulences  complexes).  Simuler  le  polluant  emis  en  suivant  les  trajectoires  de 
beaucoup  d’elements  de  fluides  «  marques  »  qui  ont  ete  emis  a  la  repartition  des  sources, 
augmente  cependant  la  difficulty  d’estimer  correctement  la  concentration  moyenne  du  polluant 
disperse,  a  partir  de  l’information  sur  la  trajectoire  de  la  particule.  La  methode  d’estimation  a 
partir  du  noyau  de  densite  a  ete  recemment  proposee  et  appliquee  avec  succes  pour  estimer  les 
concentrations  moyennes  des  modeles  de  particules  langrangiens  stochastiques.  Cependant, 

P effort  de  calcul  requis  par  cette  methode  augmente  par  un  facteur  N~  (dans  l’hypothese  ou  le 
nombre  de  lieux  recepteurs  Nr  _  auxquels  la  concentration  est  requise  est  comparable  au  nombre  de 
particules  de  fluides  Nr  utilise  dans  la  simulation  de  la  trajectoire,  tel  que  Nr  =  Np  ~  N)  et,  par 
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V 


consequent,  cette  methode,  exigeant  des  ressources  informatiques  tres  importantes,  n’a  pas  ete  tres 
utilise'e.  Nous  decrivons  ici,  un  nouvel  algorithme  qui  calcuie  P  estimation  du  noyau  du  champ  de 
la  concentration  moyenne  et  dont  la  complexity  des  calculs  est  reduite  a  N,  La  technique  utilise 
une  tessellation  (subdivision)  de  l’espace  en  cellules  cubiques  dont  la  longueur  des  cotes  h  (h  etant 
la  largeur  de  bande  de  la  fonction  du  noyau)  et  puis  associe  une  structure  de  donnees  en  liste 
chainee  avec  chaque  cellule;  cette  structure  etant  utilisee  comme  engin  de  comptabilite  suivant  les 
particules  fluides  «  marquees  »  dans  cette  cellule.  Cette  methode  rapide,  mise  au  point  ici,  a  ete 
verifiee  en  comparant  des  resultats  avec  1’ implementation  directe  de  l’estimateur  de  noyau  et  avec 
I’estimateur  classique  de  «  comptage  de  bottes  »,  pour  le  champ  de  concentration  moyenne. 

La  portee  des  resultats 

On  decrit  le  developpement  et  1’ implementation  d’un  algorithme  a  base  de  donnees  en  liste 
chainee  rapide  qui  estime  le  noyau  de  la  densite  dans  le  champ  de  concentration  moyenne,  ce  qui 
permet  a  Putilisateur  de  calculer  un  ensemble  d’estimations  du  noyau  de  la  concentration  aux  lieux 

recepteurs  Nr,  d’ordre  0(Nr)  du  travail  contrairement  a  l’ordre  0(N2r  )de  Peffort).  (Dans 
Phypothese  d’une  situation  caracteristique  ou  le  nombre  de  particules  fluides  Np  utilise  dans  la 
simulation  de  la  trajectoire  correspond  a  Np  >  Nr).  L’algorithme  rend  Putilisation  des  modeles  LS 
faisable  avec  un  ordinateur  personnel  pour  les  calculs  routiniers  de  la  dispersion,  associes  aux 
scenarios  complexes  d’emissions  d’agents  CBW,  dans  P  atmosphere. 

Les  travaux  futurs 

II  existe  plusieurs  fa9ons  de  general iser  le  code  des  algorithmes  de  donnees  a  base  de  liste  chainee 
rapide  qui  vise  a  estimer  la  densite  du  noyau,  dans  le  champ  de  concentration  moyenne,  et  /  ou  de 
rendre  ce  code  plus  efficace.  La  structure  des  donnees  a  base  de  liste  chainee  utilisee  ici  est 
appropriee  dans  le  cas  ou  la  largeur  de  bande  de  la  densite  du  noyau  associee  a  chaque  particule 
est  fixe.  En  principe,  la  largeur  de  bande  h  n’a  pas  besoin  d’etre  constante  et  dans  beaucoup 
d’applications,  il  peut  suffire  d’allouer  une  longueur  variable  souple  pouvant  dynamiquement 
s’adapter,  en  sorte  que  le  nombre  de  particules  fluides  avoisinantes,  a  l’interieur  du  support  de  la 
fonction  du  noyau  centree  sur  une  particule  donnee,  reste  constant.  La  generalisation  de 
l’estimateur  rapide  de  noyau  qui  est  propose  pour  la  largeur  de  bande  h  variable  exigera 
probablement  que  la  structure  des  donnees  a  base  de  liste  chainee  utilisee  ici  soit  remplacee  par 
une  structure  de  donnee,  un  organigramme  hierarchique,  pouvant  etre  adopte  pour  satisfaire  les 
besoins  d’une  largeur  de  bande  variable.  Dans  Papplication  de  l’estimateur  de  noyau  rapide,  le 
code  peut  etre  rendu  encore  plus  efficace  en  pre  calculant  les  valeurs  de  la  fonction  du  noyau  pour 
un  grand  nombre  de  points  representatifs  et  en  sauvegardant  les  resultats  dans  une  table  de 
recherche  de  la  valeur  de  la  fonction  pre  calculee  des  noyaux.  il  pourrait  enfin  etre  interessant 
d’explorer  la  parallelisation  du  code  de  maniere  a  ce  qu’il  puisse  etre  traite  avec  des  ordinateurs 
possedant  des  architectures  hautement  paralleles  (par  ex.  agregats  Beowulf). 


Shao.  Y.  and  Yee,  E.  (2004).  The  Rapid  Evaluation  of  Mean  Concentration  Fields  in  Lagrangian 
Stochastic  Modelling  Using  a  Density  Kernel  Estimator.  (DRDC  Suffield  TR  2004-1 86). 
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Introduction 


The  prediction  of  the  dispersion  of  contaminants  released  into  various  turbulent  flows  is  of 
practical  importance  in  a  number  of  applications  ranging  from  industrial  mixing  problems  to 
the  analysis  of  nuisance  and  hazard  in  air  quality,  pollution  and  combustion  problems 
involving  the  release  of  toxic  or  flammable  materials.  In  consequence,  much  effort  has  been 
expended  on  the  development  of  models  to  predict  mean  concentrations  of  contaminants  for  a 
given  emission-source  distribution.  However,  the  accurate  prediction  of  the  transport  and 
diffusion  of  passive  scalars  released  in  a  spatially  inhomogeneous  and  temporally 
non-stationary  turbulent  flow  (e.g.,  inhomogeneous  and  non-stationary  fields  of  mean  wind 
speed  and  turbulence  over  complex  terrain  in  the  atmospheric  boundary  layer)  is  a  challenging 
problem.  For  these  types  of  turbulent  flows,  a  Lagrangian  Stochastic  (LS)  model  describing 
the  paths  or  trajectories  of  “marked”  fluid  elements  offers  great  potential  for  studying 
contaminant  dispersion  since  they  can  model  properly  inhomogeneities,  non-steadiness,  or 
non-Gaussianity  in  the  background  velocity  distribution  (viz.,  LS  models  of  particle  motion 
are  a  particularly  attractive  option  for  the  simulation  of  dispersion  in  complex  flows). 

In  a  LS  model,  the  mean  concentration  of  a  passive  scalar  transported  by  a  turbulent  flow  field 
(whose  velocity  statistics  are  assumed  to  be  known  or  prescribed  a  priori)  is  obtained  from  the 
probability  distribution  of  the  displacements  of  independent  “marked”  fluid  elements  (or 
particles  of  tracer  which  should  not  be  viewed  as  physical  particles,  but  rather  as  markers 
keeping  track  of  the  effects  of  fluid  convection)  released  from  the  source  distribution  of  the 
scalar.  Lagrangian  approaches  for  calculation  of  dispersion  provide  a  more  natural  approach 
than  Eulerian  techniques  (usually  based  on  eddy-diffusivity  concepts  which  are  strictly  valid 
only  for  the  case  where  the  characteristic  length  scale  of  the  spatial  gradients  in  the  mean  of  a 
scalar  quantity  is  larger  than  the  integral  length  scale  of  the  turbulence,  implying  that  these 
techniques  cannot  be  applied  to  the  simulation  of  dispersion  close  to  a  source).  Moreover,  LS 
models  can  be  designed  rigorously,  based  on  a  small  set  of  principles,  to  give  reliable 
simulations  of  dispersion  in  complex  inhomogeneous,  non-stationary  or  non-Gaussian  flows 
where  other  methods  based  on  eddy  diffusivity  concepts  and/or  other  semi-empirical 
approaches  (e.g.,  similarity  theory,  asymptotic  solutions,  etc.)  are  invalid  or  inadequate.  In  this 
regard,  the  well-mixed  condition  due  to  Thomson  [1]  provides  the  most  rigorously  correct 
theoretical  criterion  for  the  formulation  of  LS  models  of  turbulent  dispersion. 

In  a  stochastic  or  random-walk  particle  model,  dispersion  is  simulated  by  following  the 
trajectories  of  a  large  number  of  “marked”  fluid  particles  released  from  the  source  into  the  flow 
domain,  with  the  mean  concentration  of  the  scalar  determined  from  the  statistics  of  the  particle 
displacements.  This  gives  rise  to  a  primary  difficulty  in  the  application  of  LS  models;  namely, 
the  “correct”  and  numerically  efficient  estimation  of  the  mean  concentration  at  a  given  location 
and  time  from  the  information  embodied  in  the  “marked”  particle  trajectories.  A  stochastic  or 
random  walk  model  of  particle  motions  is  a  grid-free  method  in  which  particles  evolve  in  time 
according  to  stochastic  equations  of  motion.  Nevertheless,  it  is  common  practice  in 
atmospheric  dispersion  modelling  with  LS  models  to  calculate  the  mean  concentration  on  a 
fixed  grid  using  a  particle-in-cell  or  box-  counting  method  [2,  3, 4,  5],  In  the  box-counting 
method,  the  mean  concentration  is  proportional  to  the  local  particle  number  density  (since  the 
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mass  carried  by  each  particle  is  unchanged  along  the  trajectory  of  the  particle).  More 
specifically,  the  mean  concentration  is  estimated  by  multiplying  the  number  of  particles  in  a 
grid  box  by  the  mass  carried  by  each  particle  (assumed  equal),  and  dividing  this  quantity  by 
the  volume  of  the  grid  box.  Unfortunately,  the  box-counting  method  for  mean  concentration 
estimation  depends  on  the  choice  of  the  size  and  position  of  the  grid  boxes  used  for  the 
estimation.  To  overcome  this  problem,  Yamada  et  al.  [6],  Uliasz  [7],  and  de  Haan  [8]  used  a 
density  kernel  method  for  mean  concentration  estimation  which  fits  well  into  a  LS  particle 
model  framework  since  the  mean  concentration  is  calculated  without  reference  to  a  grid.  A  LS 
particle  model  coupled  with  the  use  of  a  density  kernel  estimator  allows  the  mean 
concentration  to  be  determined  without  reference  to  a  grid  (and,  hence,  provides  a  truly 
grid-free  Lagrangian  Monte  Carlo  method). 

In  spite  of  the  advantages  of  the  kernel  density  estimator  and  the  deficiencies  in  the 
box-counting  method,  the  latter  method  is  still  predominantly  used  in  atmospheric  dispersion 
modelling  for  the  estimation  of  the  mean  concentration  field.  This  is  because  the  kernel 
density  estimator  for  the  mean  concentration  is  computationally  much  more  demanding  than 
the  simpler  box-counting  method.  More  specifically,  when  implemented  in  a  straightforward 
manner  the  computational  work  of  the  kernel  density  estimator  scales  as  0(NpNr)  (Np  and  Nr 
being  the  number  of  particles  used  in  the  simulation  and  the  number  of  receptor  points  at 
which  the  mean  concentration  is  calculated,  respectively),  whereas  in  the  box-counting  method 
the  computational  work  is  of  order  0(Nr).  Since  large  numbers  of  particles  are  typically 
required  to  accurately  model  the  mean  concentration  field,  the  overall  computational  cost  of 
using  the  kernel  density  estimator  in  a  LS  particle  model  would  quickly  become  prohibitive.  In 
consequence,  it  would  be  desirable  to  develop  an  algorithm  for  mean  concentration  estimation 
using  density  kernels  for  which  the  computational  work  scales  as  0{Nr),  enabling  a 
computationally  efficient  and  feasible  kernel  density  method  for  estimation  of  concentrations 
from  a  stochastic  Lagrangian  particle  model.  It  is  the  purpose  of  this  report  to  develop  and 
implement  an  algorithm  which  can  calculate  the  mean  concentration  from  a  disordered  set  of 
particle  positions  (obtained  from  a  LS  particle  model)  in  0{Nr)  computational  work  as 
opposed  to  0(NpNr)  work.  The  reduction  of  the  computational  complexity  using  our  new 
algorithm  will  allow  kernel  density  estimators  for  mean  concentration  to  be  used  routinely  in 
LS  models  for  atmospheric  dispersion. 


Mean  Concentration  Estimation 


The  LS  modelling  of  the  dispersion  of  a  passive  tracer  consists  of  simulating  numerically  the 
motion  of  many  particles  of  tracer  (“marked”  fluid  elements)  released  from  the  source 
distribution,  in  order  to  build  up  a  picture  of  the  concentration  distribution.  More  specifically, 
if  S(x',t')  denotes  the  spatial-temporal  density  distribution  of  the  source  [viz.,  the  mass  of 
source  material  per  unit  volume  per  unit  time  interval  at  the  source  space-time  point  (x',f')3> 
then  the  mean  concentration  C(x,t)  at  the  receptor  space-time  point  (x,t)  is  determined  from 

[9] 

(1)  C(x,f)  =  f  [  pLfatlx'^Six'j^dt'dx', 

JrJ  o 
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where  pi(x,t\x',t')  =  (8 (x  — X(/;x/,/')))  is  the  conditional  probability  density  function 
(PDF)  that  a  marked  particle  released  at  the  source  space-time  point  (x',t')  will  be  found  at 
receptor  (observation)  space-time  point  (x,t).  Here,  S(-)  is  the  Dirac  delta  function,  (•)  is  used 
to  denote  the  averaging  over  the  samples  of  the  turbulent  velocity  field,  and  X(t;x',t')  denotes 
the  Lagrangian  spatial  coordinates  defining  the  trajectory  in  physical  space  of  a  particle 
“marked”  at  the  source  point  (x',t').  In  Eq.  (1),  the  spatial  integration  is  over  the  domain  (or, 
volume)  V  that  contains  the  emission-source  distribution  S(x',t'). 

In  accordance  to  Eq.  (1),  to  estimate  the  mean  concentration  at  receptor  point  x  and  time  t 
requires  the  representation  of  the  integral  here  as  an  expectation  of  a  random  estimator  defined 
on  Lagrangian  trajectories.  Since  the  exact  form  of  pi(x,/|x'.r')  is  unknown,  we  need  to 
approximate  this  particle  displacement  PDF  from  the  solutions  of  stochastic  evolution 
equations  describing  trajectories  of  fluid  particles  in  a  turbulent  flow.  The  evolution  equations 
for  the  position  X(r)  and  velocity  U(t)  of  a  tagged  fluid  particle  in  a  Lagrangian  Stochastic 
model  is  typically  based  on  a  continuous  Markovian  evolution  in  position-velocity  (phase) 
space  specified  by  a  stochastic  differential  equation  of  the  form  [1] 


(2) 

dUi  =  a,{V,X,t)dt  +  (C0e)1/24/(0 

and 

(3) 

dXj  =  Ufdt, 

where  dU,  and  dX{  ( i  —  1,2,3)  denote  the  change  (or,  increment)  in  velocity  and  position  over 
the  time  interval  dt,  respectively.  More  specifically,  dUt  =  Ut{t  +  dt)  —  Ut{t )  is  the 
infinitesimal  increment  of  the  velocity  Ut  following  the  marked  fluid  particle.  Appearing  in  the 
equations  is  the  drift  coefficient  vector  u,(U,X,t),  the  universal  constant  C0  associated  with  the 
Lagrangian  velocity  structure  function,  the  mean  rate  of  dissipation  of  turbulence  kinetic 
energy  e,  and  the  isotropic  vector  incremental  Weiner  process  d^j(t).  The  value  Co  =  3.5  is 
recommended  for  the  LS  model  used  here  [2,  10].  The  increments  d^  have  a  joint-normal 
distribution  with  zero  means  and  an  isotropic  covariance  matrix, 

(4)  (dUt))  =  0,  (dUt)d^j(t  +  x)>  =  dt5ijb(x), 

where  8 y  denotes  the  Kxonecker  delta  function.  Finally,  it  should  be  emphasized  that  in  a 
“correctly”  formulated  LS  model,  the  drift  coefficient  vector  should  be  chosen  to  satisfy  the 
well-mixed  condition  [1]  which  guarantees  that  the  model  will  maintain  the  correct 
phase-space  distribution  of  particles  in  the  sense  that  once  the  particles  of  tracer  are 
well-mixed  in  the  flow,  they  will  remain  so. 

The  conventional  method  for  evaluating  the  mean  concentration  of  Eq.  (1)  is  to  apply  the 
box-counting  method  and  count  the  number  of  particles  in  a  discrete  grid  cell  centered  at  x  at 
time  t.  This  process  is  identical  to  the  calculation  of  a  three-dimensional  histogram.  In 
particular,  the  estimate  C  for  the  mean  concentration  at  a  given  location  x  and  time  t  is  given  by 

(5)  C(x,t)  —  — x  Am. 

A/ 

where  Npi(t)  is  the  number  of  “marked”  fluid  particles  at  time  t  in  the  spatial  cell  /  (that 
contains  the  receptor  position  x),  A /  is  the  volume  of  cell  /,  and  Am  is  the  mass  of  tracer  carried 
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by  each  fluid  particle.  Each  fluid  particle  in  the  ensemble  of  Np  particles  is  assumed  to  carry  a 
constant  mass  Am.  Notice  that,  for  a  given  spatial  cell  /,  the  sample  size  is  usually  relatively 
small  (i.e.,  Npi  <C  Np),  and  hence  the  statistical  error  in  the  concentration  estimate  C(x.t)  will 
generally  be  relatively  large  with  the  error  scaling  as  Np  .In  consequence,  box-counting 
estimates  for  the  mean  concentration  will  therefore  generally  require  large  sample  sizes  (i.e., 
Np  large)  if  the  statistical  error  is  to  remain  small.  Alternatively,  it  is  possible  to  reduce  the 
statistical  error  by  increasing  the  size  (volume)  of  the  spatial  cell,  but  this  may  result  in  the 
structural  details  in  the  mean  concentration  being  oversmoothed  (leading  to  a  larger  bias  in  the 
mean  concentration  estimate).  In  summary,  many  particles  per  spatial  cell  are  required  to  keep 
the  bias  and  statistical  error  in  the  mean  concentration  estimate  small,  but  this  requirement 
makes  trajectory  simulations  computationally  expensive  and  hence  prohibitive. 

For  these  reasons,  Uliasz  [7]  and  de  Haan  [8]  proposed  another  method  for  the  estimation  of 
the  mean  concentration  based  on  density  kernels.  The  foundation  of  the  kernel  estimate  of  a 
field  variable  at  a  point  is  interpolation  theory  generalized  to  the  case  of  interpolation  from  the 
disordered  positions  of  an  ensemble  of  marked  particles,  where  the  underlying  turbulent 
velocity  field  “deals”  out  the  disordered  positions.  Within  the  volume  of  the  physical  space  V, 
the  mean  concentration  of  the  dispersing  scalar  is  represented  by  Np  particles,  each 
representing  a  mass  of  tracer  A m  with  the  n-th  particle — the  numbering  being 
arbitrary — having  a  position  X^(t)  and  velocity  IjM(t)  at  time  t.  Then  the  mass  density 
function  (concentration)  C*(x,t )  (which  can  be  interpreted  to  be  the  discrete  representation  of 
the  mean  concentration  distribution  C(x,t)  of  a  dispersing  scalar  and,  as  such  can  be 
considered  to  be  an  estimate  of  this  distribution)  is  defined  by 

Np 

(6)  C*(x,0  =  Am^5(x-xW(0), 

71=1 

where 

(7)  5(x-xW(/))  =  6(x,  -4n)(t))8(x2 -X,('°(0)5(x3 

is  the  3-dimensional  Dirac  delta  function.  For  simplicity  (but  without  any  loss  in  generality), 
assume  that  the  tracer  is  released  from  an  instantaneous  point  source  located  at  the  origin 
x  =  0  with  the  release  occurring  at  time  t  =  0.  Denoting  by  Q  the  total  mass  of  tracer  released 
(Q  =  Np&m),  the  (mathematical)  expectation  of  Eq.  (6)  is 

(8)  <C*(x,t)>  =  0(8(x-X*(O)}, 

where  *  refers  to  any  particle  n  (1  <  n  <  Np).  Comparing  Eq.  (8)  with  Eq.  (1)  for  the  case  of 
an  instantaneous  point  source  release  at  x  =  0  and  t  =  0  {Six' ,t')  =  05(x/)5(t/)),  it  is 
straightforward  to  show  that  the  mass  density  function  in  Eq.  (6)  is  a  consistent  (unbiased) 
representation  of  the  mean  concentration  in  the  sense  that 

(9)  <C*(x,f))=C(x,f). 


4 
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The  use  of  an  interpolation  function  that  provides  the  kernel  estimate  of  the  mean 
concentration  at  a  space-time  point  involves  replacing  the  3-dimensional  Dirac  delta  function 
in  Eq.  (6)  by  a  weighting  function  (or,  “kernel”)  that  defines  how  much  of  the  mass  carried  by 
a  “marked”  fluid  particle  contributes  to  the  mean  concentration  field  at  the  point  x  at  time  t. 
The  kernel  estimate  of  C(x.t)  is  then  given  by 

Np  1 

(10)  c(x, o  =  Am  ]T-x(x-xM  (;);/*), 

n=  1  n 

where  h  is  the  width  of  the  (interpolating)  kernel  function  K(r,h).  The  kernel  function  is  a 
non-negative  valued  function  [ K(r;h )  >  0]  that  satisfies  the  property 

(11)  [  K(v,h)dr=l, 

J  a.s. 

where  the  integral  is  taken  over  all  space  (a.s.). 


The  kernel  method  does  not  require  the  physical  space  to  be  partitioned  into  grid  cells  and 
produces  a  smooth  mean  concentration  field  with  a  much  smaller  number  of  particles  than  is 
required  for  the  box-counting  method.  The  width  (or,  smoothing  length)  h  determines  the 
degree  of  smoothing  of  the  mean  concentration  field.  The  smaller  values  of  h  give  a  more  local 
estimate  for  the  mean  concentration,  but  result  in  fewer  particles  giving  significant 
contributions  at  a  given  (receptor)  point  x  and,  hence,  more  statistical  error  in  the  estimate. 
Taking  h  larger  allows  more  “marked”  fluid  particles  to  contribute  to  the  mean  concentration 
estimate  at  a  receptor  location  x  reducing  the  statistical  error,  but  taking  h  too  large  will  smooth 
out  detailed  features  from  the  mean  concentration  distribution  increasing  the  bias  in  the 
estimate.  In  consequence,  the  choice  of  the  bandwidth  h  is  a  compromise  between  smoothing 
enough  and  not  smoothing  too  much  to  smear  out  the  real  features  in  the  mean  concentration 
distribution.  Mathematically,  there  is  a  compromise  between  the  bias  and  variance  of  C(x,t), 
which  increases  and  decreases,  respectively,  as  h  is  increased.  Theory  for  density  kernel 
estimators  [11]  suggests  that  the  kernel  bandwidth  should  be  chosen  proportional  to  Npl  5 ,  but 
the  constant  of  proportionality  depends  on  the  unknown  mean  concentration  distribution. 


Following  de  Hann  [8],  we  consider  two  types  of  kernel  functions  K  with  finite  (or,  compact) 
support.  The  latter  implies  that  only  a  subset  of  the  “marked”  fluid  particles  contributes  to  each 
kernel  estimate  of  the  mean  concentration.  In  particular,  we  use  kernel  functions  with  the 
general  form 


(12) 


Q,a(l  —S2y  ; 

0; 


s2  <  1, 
s2  >  1, 


1  U 

where  s  =  (r  •  r)  '  /h,  a  is  an  exponent,  d  is  the  number  of  dimensions,  and  Q.a  is  the 
normalization  constant.  In  this  report,  we  consider  kernel  functions  with  a  =  1  (parabolic 
kernel,  which  is  referred  to  also  as  Epanechnikov  kernel  [12])  and  a  =  4  (quadweight  kernel). 
For  these  two  kernels,  the  normalization  constant  Q.a  [obtained  from  the  normalization 
constraint  Eq.  (11)]  takes  the  values  £  or  for  a  —  1  and  |  or  for  a  —  4  in  two-  and 
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three-dimensions  (d  =  2  and  3),  respectively.  One  of  the  key  advantages  of  using  Eq.  (12),  as 
opposed  to,  say  a  Gaussian  kernel  function  (K(r;h)  =  exp(—s2/2) / (2nh)d/2),  is  that  it  has  a 
finite  support  so  that  “marked”  fluid  particle  contributions  to  mean  concentration  at  a  given 
(receptor)  point  x  are  exactly  zero  for  particles  at  distances  r  =  (r  ■  r)1/2  >  h  (or,  equivalently, 

5  =  rjh  >  1).  Figure  1  shows  the  profiles  of  non-normalized  parabolic  and  quadweight  kernels 
(i.e.,  profiles  of  K(r,h)/Cd.a  =  K(s')). 

A  Fast  Algorithm  for  Kernel  Density  Estimation _ 

Algorithm  description 

The  naive  implementation  of  the  kernel  estimate  [cf.  Eq.  (10)]  for  the  mean  concentration  field 
is  simple.  Since  the  kernel  functions  considered  here  have  a  compact  support  [refer  to 
Eq.  (12)],  only  a  subset  of  the  fluid  particles  will  contribute  to  the  kernel  estimate  at  a  given 
(fixed)  receptor  location  x,  corresponding  as  such  to  the  nearest  neighboring  particles  for  the 
receptor  position.  The  naive  evaluation  of  Eq.  (10)  then  consists  of  an  all-pair  search  approach 
whereby  the  distance  =  ((x-XW(t))  •  (x-XM(r)))1^2  ( n  =  1,2  ,...,NP)  of  each  tagged 
fluid  particle  from  the  receptor  location  x  is  calculated,  and  if  this  distance  is  smaller  than  h 
(bandwidth  of  the  kernel  function  centered  about  the  particle  position,  which  also  defines  the 
support  domain  of  the  function),  the  mass  carried  by  the  particle  contributes  to  the  mean 
concentration  at  receptor  point  x.  The  all-pair  search  approach  needs  to  be  carried  out  for  each 
receptor  location  xr  (r—  \  ,2,...,Nr,  where  Nr  is  the  number  of  receptor  locations  where  the 
mean  concentration  is  calculated),  and  at  each  of  these  locations  the  searching  is  performed  for 
all  Np  particles  at  the  given  time  t.  It  is  clear  that  the  complexity  of  the  all-pair  search 
algorithm  for  the  computation  of  the  kernel  estimate  of  the  mean  concentration  field  is  of  order 
0{NrNp)  0(N2)  (if  Nr&Np  =  N).  Furthermore,  this  naive  algorithm  needs  to  be  applied  at 
all  time  steps  (since  the  tagged  fluid  particle  positions  are  evolving  in  time)  so  that  the 
computational  cost  of  the  conventional  (albeit  straightforward)  implementation  of  the  kernel 
estimate  of  the  mean  concentration  field  is  computationally  prohibitive  (and,  typically,  too 
expensive  to  be  applied  in  practical  problems  where  the  number  of  particles  is  large). 

The  straightforward  implementation  of  the  kernel  estimate  for  the  mean  concentration 
corresponds  to  a  huge  waste  of  computational  time  since  only  a  small  fraction  of  the  total  Np 
fluid  particles  will  contribute  to  the  concentration  at  a  given  (fixed)  receptor  location.  Our 
approach  for  improving  the  computational  efficiency  of  the  kernel  estimator  (grid-free  method) 
is  to  maintain  a  neighbor  list  for  each  fluid  particle.  The  methodology  used  to  keep  track  of 
nearest  neighbor  fluid  particles  relative  to  a  given  receptor  point  is  to  divide  the  physical 
volume  V  (flow  domain)  into  cells  or  elements  (viz.,  a  temporary  mesh  is  overlaid  on  the  flow 
domain  to  provide  a  tesselation  of  the  domain).  Each  cell  in  the  tesselation  of  the 
3-dimensional  domain  is  a  cube  (or,  a  square  for  a  2-dimensional  domain)  with  side  length 
equal  to  h  (smoothing  length  of  the  kernel  function),  the  reason  being  that  the  kernel  functions 
centered  on  the  fluid  particles  only  contribute  to  the  concentration  at  a  receptor  location  if  the 
distance  r  between  the  particle  position  and  the  receptor  location  is  less  than  h  (i.e,,  r  <  h).  In 
consequence,  with  cells  of  side  length  h,  we  only  need  to  check  the  current  cell  which  contains 
the  receptor  location  and  nearest  neighboring  cells  in  order  to  find  all  possible  fluid  particles 
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that  can  contribute  to  the  mean  concentration  at  the  receptor  location. 

To  facilitate  using  cells  as  a  book-keeping  device  to  keep  track  of  fluid  particle  positions,  each 
cell  in  the  superimposed  structured  grid  is  associated  with  a  linked  list.  Each  “marked”  fluid 
particle  in  the  flow  domain  is  assigned  to  the  cell  it  resides  in  and  identified  through  the  linked 
list  associated  with  that  cell.  The  cells  in  the  tesselation  of  the  physical  domain  can  be 
identified  by  the  index  set  (i,j)  or  ( i,j,k ),  respectively,  for  a  2-dimensional  or  3-dimensional 
domain.  If  ( i,j,k )  is  the  index  set  for  cell  P  in  a  3-dimensional  physical  region,  then  the  cells 
immediately  to  the  right  and  left  of  P  are  labelled  (i  +  1  ,j,  k)  and  (i  —  1  ,j,  k),  the  cells 
immediately  in  front  and  back  are  (i,j+  1  ,k)  and  (i,j  —  l,fc),  and  the  cells  directly  above  and 
below  are  (i,  j,k 4- 1 )  and  (i,j,k  —  1).  For  a  given  receptor  location  x  contained  in  cell  P  with 
label  (i,j, k ),  all  fluid  particles  which  can  contribute  to  the  concentration  at  x  can  only  reside  in 
cell  P  and/or  in  the  one  of  its  immediately  adjoining  (nearest  neighbor)  cells  with  labels 
(i±a,ydhP,fc±y),  where  a,  (3,  y€  {0,1}  excluding  the  choice  a  =  p  =  y=  0  (which 
corresponds  to  cell  P  itself).  Hence,  for  the  computational  cell  (box)  P  labelled  by  the  index 
(i i,j,k ),  the  list  of  all  nearest  neighbor  cells  to  P,  denoted  by  near (P),  is  the  set  of  all  cells 
labelled  by  the  index  set  {(i±a,j  ±p,£±y)|a,  P,  y  G  {0,  l},a  p  ^  y  =  0}.  Therefore,  the 
search  for  fluid  particles  that  can  contribute  to  the  concentration  at  x  is  confined  to  only  9  or  27 
cells  in  2-  or  3-dimensional  space,  respectively.  The  linked  list  associated  with  each  cell  allows 
the  fluid  particles  assigned  to  each  cell  to  be  chained  together  for  easy  access. 

The  superimposed  grid  allows  the  “binning”  (or  sorting)  of  the  fluid  particles  into  Nc  cells, 
with  the  average  number  of  particles  per  cell  being  about  Np/Nc  (assuming  the  particles  are 
uniformly  distributed  in  the  physical  domain).  If  the  average  number  of  fluid  particles  per  cell 
is  small  (Np  »  Np/Nc,  the  latter  of  which  is  of  order  0(1)),  the  computational  complexity  of 
the  linked  list  algorithm  for  the  kernel  estimate  of  the  mean  concentration  field  is  of  0{Nr). 
This  computational  effort  is  optimal  since  the  program  must  at  least  loop  through  all  Nr 
receptor  locations  in  order  to  compute  the  mean  concentration  field  at  a  particular  time  t  (viz., 
traversing  a  “data”  set  of  Nr  receptor  locations  has  a  minimal  computational  complexity  of 
order  0(Nr)). 


Programming  considerations 

A  computer  program  that  implements  the  linked  list  (fast)  algorithm  for  the  kernel  estimation 
of  the  mean  concentration  field  has  been  written.  The  program  contains  less  than  several 
hundred  lines  of  Fortran  90/95  code.  In  this  subsection,  we  summarize  some  of  the  more 
technical  details  of  our  implementation  of  the  fast  kernel  density  estimator. 

The  beginning  of  the  code  shown  below  contains  definitions  of  the  components  for  one  derived 
data  type:  namely,  cells.  This  derived  data  type  cells  contains  an  integer  variable  counter 
to  hold  the  fluid  particle  number  index,  four  real  variables  mass,  x,  y,  and  z  to  hold  the  fluid 
particle  mass  and  the  coordinates  (x,y,z)  of  the  fluid  particle  position  in  3-dimensional  space, 
and  a  pointer  variable  NextParticle  to  a  derived  data  type  cells  variable  as  its  last 
component.  This  pointer  points  to  the  next  item  (particle)  in  the  linked  list  of  particle 
properties  contained  in  a  given  computational  cell.  A  linked  list  of  particles  associated  with  a 
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given  cell  is  simply  a  series  of  variables  in  the  derived  data  type  cells  with  the  pointer  from 
each  derived  data  type  variable  pointing  to  the  next  variable  in  the  list.  The  pointer 
Next  Particle  of  the  last  element  of  the  linked  list  is  set  to  NULL,  a  convenient  sentinel 
marking  the  end  of  the  list. 

Each  cell  of  the  tesselation  of  the  physical  region  has  a  linked  list  associated  with  it,  the  latter 
used  simply  as  a  convenient  book-keeping  device  to  keep  track  of  the  particles  inside  that  cell 
at  a  given  time  t.  This  feature  is  very  useful  because  it  permits  us  to  construct  various  types  of 
dynamic  data  structures  linked  together  by  successive  pointers  during  the  execution  of  a 
program.  We  use  a  linked  list  to  hold  the  fluid  particle  properties  in  a  cell  because  the  number 
of  fluid  particles  in  a  given  cell  can  vary  with  time  as  the  position  and  velocity  of  the  particles 
evolve  according  to  the  stochastic  differential  equations  exhibited  in  Eqs.  (2)  and  (3),  and  the 
linked  list  permits  us  to  add  (remove)  elements  to  (from)  the  list  one  at  a  time  without  knowing 
in  advance  how  many  elements  will  ultimately  be  in  the  list.  Since  there  are  mutiple  cells  in 
the  flow  domain,  another  derived  data  type  cell-pointer  is  defined  and  used  in  the 
declaration  of  3-dimensional  arrays  of  pointers  to  derived  data  type  cells.  Here, 
cell-pointer  is  defined  to  be  a  pointer  to  the  derived  data  type  cells.  Finally,  two 
3-dimensional  arrays  of  pointers  (head  and  tail)  to  variables  of  derived  data  type  cells  are 
also  defined  to  point  to  the  first  and  last  variables  in  the  linked  list  for  each  cell.  To  access  an 
element  in  the  linked  list,  we  start  at  the  head  of  the  list  and  use  the  pointers  NextParticle  of 
successive  elements  to  move  from  element  to  element  until  the  desired  element  is  reached. 
With  the  singly-linked  lists  used  here,  the  list  can  be  traversed  in  only  one  direction  (from  head 
to  tail)  because  each  element  contains  no  link  to  its  predecessors. 

The  following  is  a  formal  definition  of  some  derived  data  types  required  by  the  algorithm. 

TYPE  : :  cells 

INTEGER  : :  counter 
REAL  : :  mass 
REAL  : :  X 
REAL  ::  y 
REAL  : :  Z 

TYPE  (cells),  POINTER  ::  NextParticle 

END  TYPE 

TYPE  celljpointer 

TYPE (cells),  POINTER  ::  pptr 

END  TYPE  celljpointer 

TYPE (celljpointer) ,  ALLOCATABLE,  DIMENSION^,:,:)  ::  head 

TYPE (celljpointer) ,  ALLOCATABLE,  DIMENSION (:,:,: )  ::  tail 

TYPE  (cells)  ::  temp  !  Temporary  variable 

Having  defined  some  useful  data  structures,  the  fast  algorithm  for  kernel  estimation  of  the 
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mean  concentration  field  consists  of  the  following  main  steps: 


1 .  Tesselate  the  flow  domain  into  cells  (cubes)  with  side  length  h  (bandwidth  of  the  kernel 
function). 

2.  Open  the  input  file  containing  the  “marked”  fluid  particle  properties  (index,  mass, 
location)  at  a  given  time  t. 

3.  Read  in  the  properties  (index,  mass,  location)  of  each  fluid  particle,  determine  the  cell  the 
particle  resides  in,  and  insert  the  relevant  particle  information  into  the  linked  list  for  that 
cell. 

4.  Define  the  set  of  receptor  locations  at  which  to  calculate  the  mean  concentration. 

5.  For  each  receptor  location,  determine  the  cell  P  that  it  resides  in,  extract  the  fluid  particle 
properties  from  the  linked  list  for  this  cell,  and  calculate  the  contribution  of  each  of  these 
fluid  particles  to  the  kernel  density  estimator  of  the  mean  concentration. 

6.  Identify  the  26  nearest  neighbor  cells  (near(P))  to  the  cell  P  containing  the  receptor 
location,  extract  the  fluid  particle  properties  from  the  linked  lists  for  each  of  these 
adjoining  cells,  and  for  each  particle  calculate  the  Euclidean  distance  r  between  the 
particle  and  receptor  location.  For  each  particle  whose  distance  r  from  the  receptor 
position  is  less  than  h,  calculate  the  contribution  of  this  particle  to  the  kernel  density 
estimator  of  the  mean  concentration. 

7.  Output  the  mean  concentration  field  at  the  given  time  t. 

Following  is  a  formal  description  of  the  algorithm. 


Algorithm 


Input  functional  form  of  kernel  function 
Input  bandwidth  h  of  kernel  function 

Subdivide  the  flow  domain  into  contiguous  cells  of  side  length  h 
Open  file  containing  fluid  particle  property  information 
WHILE 

Read  fluid  particle  property  information  into  temp 
IF  read  not  successful  EXIT 

Determine  index  of  cell  where  current  particle  resides 
Get  linked  list  for  current  cell 
!  determine  which  cell  ( ip , j p , kp)  the  particle  is  located  in 

ip  =  INT((temp%x  -  Xmin) /bandwidth)  +  1 


DRDC  Suffield  TR  2004-186 


9 


jp  =  INT((temp%y  -  Ymin) /bandwidth)  +  1 
kp  =  INT((temp%z  -  Zmin) /bandwidth)  +  1 

IF  head(ip, jp,kp)  (of  linked  list)  is  not  associated  THEN 

!  The  list  is  empty 

ALLOCATE  head(ip, jp,kp) 

!  Tail  points  to  first  value 

tail (ip, jp,kp)  =>  head(ip,jp,kp) 

I  Nullify  pointer  NextParticle  within  1st  value 
!  since  there  is  nothing  to  point  to  yet 

NULLIFY  tail  (ip, jp,kp) INextParticle 

!  Store  particle  information  (index, mass, location) 

!  in  linked  list 

tail (ip, jp,kp)  <--  temp 

ELSE 

!  The  list  already  has  particle  information  inserted 
ALLOCATE  tail (ip, jp, kp) %NextParticle 
!  Tail  now  points  to  new  last  value 

tail (ip, jp,kp)  =>  tail (ip, jp,kp) %NextParticle 
I  Nullify  pointer  within  new  last  value 
NULLIFY  tail  (ip, jp,kp) %Nextparticle 
!  Store  new  particle  information 
tail (ip, jp,kp)  <--  temp 
END  of  IF 
END  of  WHILE 

Define  receptor  locations  and  calculate  index  (ic,jc,kc) 
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of  cell  where  each  receptor  resides 


DO  for  each  receptor  location  with  index  (ic,jc,kc) 

Get  linked  list  for  receptor  cell 

ptr  =>  headfic, jc,kc)  !  go  back  to  the  head  of  the  list 
WHILE  ptr  is  associated 

Get  the  particle  coordinates 

Calculate  contribution  to  concentration  at  receptor 
ptr  =>  ptr%Nextpointer 

END  of  WHILE 

Get  linked  lists  for  all  the  nearest  neighbor  cells  to 
receptor  cell  P  with  index  (ic,jc,kc) 

DO  for  each  linked  list  in  near(P) 

ptr  =>  head (near (P) )  !  go  back  to  the  head  of  the  list 

WHILE  ptr  is  associated 

Get  the  particle  coordinates 

Calculate  Euclidean  distance  from  particle  to  receptor 
If  distance  less  than  h,  calculate  contribution 
to  concentration  at  receptor 
ptr  =>  ptr%Nextpointer 
END  of  WHILE 
END  of  DO 

END  of  DO  !  next  receptor  location 

WRITE  out  mean  concentration  field  at  all  receptor  points 


A  Fortran  90/95  computer  program  that  implements  the  fast  kernel  estimator  for  the  mean 
concentration  is  available  from  the  authors  upon  request.  The  computer  program  has  been 
implemented  using  both  2-dimensional  and  3-dimensional  kernels. 

Simple  Lagrangian  Stochastic  Particle  Model 


As  a  test  of  our  new  method  for  the  fast  kernel  estimation  of  the  mean  concentration,  we  have 
written  a  simple  trajectory-simulation  model  for  turbulent  dispersion  from  an  instantaneous 
and  continuous  point  source  in  a  neutrally  stratified  inhomogeneous  shear  flow  (constant  stress 
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layer  in  the  atmospheric  boundary  layer).  The  trajectory-simulation  approach  to  modelling 
dispersion  of  a  passive  tracer  consists  of  simulating  the  trajectories  of  many  particles  of  tracer 
to  build  up  a  picture  of  the  mean  concentration  distribution.  In  this  section,  we  use  Cartesian 
tensor  notation,  but  for  convenience  ocassionally  we  also  use  Cartesian  coordinates  (x.y.z)  and 
(m,v,w)  to  denote  the  streamwise,  cross-stream,  and  vertical  coordinates  and  velocities, 
respectively. 


Particle  model  description 

The  dispersion  of  a  passive  scalar  is  described  by  the  stochastic  differential  equation 
(generalized  Langevin  equation)  exhibited  in  Eqs.  (2)  and  (3).  The  drift  term  (or,  conditional 
particle  acceleration)  a,-  in  Eq.  (2)  is  constrained  for  a  particular  flow  by  specifying  the  the 
Eulerian  velocity  statistics  for  that  flow,  and  requiring  that  the  Eulerian  statistics  determined 
from  the  LS  model  [Eqs.  (2)  and  (3)]  reproduce  the  specfied  statistics.  A  sufficient  condition 
for  this  to  be  achieved  is  embodied  in  the  well-mixed  criterion  (Thomson  [1])  which  requires 
the  Eulerian  velocity  PDF  PE(u;x,t)  satisfy  the  Fokker-Planck  equation  for  the  joint  PDF  of 
velocity  and  position  P(u,x,t)  associated  with  the  stochastic  differential  equation  system 
specfied  by  Eqs.  (2)  and  (3).  Application  of  the  well-mixed  criterion  constrains  a,-  as 

1  dP 

(13)  aiPE  =  -C0£ +  <>,, 

where 

3<t>/  _  dPE  dujPE 
U  j  du,  ~  dt  dx,  ’ 

and  <j>,  — >  0  as  |u|  — *■  °°  (|  •  |  denotes  the  Euclidean  norm).  Note  that  Eq.  (14)  does  not  determine 
<{>,  uniquely  since  an  arbitrary  solenoidal  vector  function  in  velocity  space  (i.e.,  function  whose 
divergence  with  respect  to  u,  vanishes)  can  be  added  to  <j),  and  still  satisfy  Eq.  (14).  Hence,  the 
well-mixed  criterion  allows  one  to  determine  ‘a’  drift  function  rather  than  ‘the’  drift  function 
since  one  can  add  any  solenoidal  vector  function  \| /,■  (<A| /,  /3w;  =  0)  to  cj),  without  changing 
dcf), /diii  and  the  interpretation  of  this  “gauge”  freedom  in  terms  of  the  statistical  properties  of 
the  flow  is  awkward. 


We  consider  the  case  of  Gaussian  stationary,  inhomogeneous  turbulence  where  PE  is  given  by 
(15)  PE(u;x,t)  =  ^jexp  (~(«,-£7/)r0-(u7-^  , 


where  T;y  =  is  the  inverse  of  the  Reynolds  stress  tensor  Vy  (i.e.,  Vy  =  («•«'•)  is  the 

covariance  between  the  j-th  and  y'-th  components  of  the  Eulerian  velocity  where  a  primed 
quantity  refers  to  the  deviation  from  the  mean  value  for  that  quantity),  y  is  the  determinant  of 
T/y,  and  £/,  is  the  mean  velocity.  Summation  over  repeated  tensor  indices  i,j,k, ...  is  implied 
(Einstein  summation  convention).  In  this  case,  Thomson  [1]  provided  the  following  form  for  a 
conditional  acceleration  vector  <j>,-  (with  u\  =  —  Uj  being  the  velocity  fluctuation)  that  is 

consistent  with  the  well-mixed  criterion: 


(16) 


_  fj  Mi  dUt  .ldVa  l  -  dVn  1  dVtl  ,  , 
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The  expanded  form  of  Eq.  (16)  is  complicated.  In  this  report,  we  consider  a  simple 
neutrally-stratified  flow  with  mean  velocity  vector  U  =  ((7,0,0)  and  turbulence  velocity 
statistics  varying  only  in  the  z-direction;( vertical  direction).  Furthermore,  we  make  the 
assumption  that  we  can  ignore  the  covariance  between  the  different  velocity  components  (i.e., 
{u'ju'j)  =  0  for  i  /  j).  In  this  case,  Thomson’s  solution  for  the  drift  coefficient  reduces  to  with 
i  =  1,2,3  =  u,v,w: 


(17) 


au  = 


av  = 


aw  = 


Coe.  ,  ,  dU  ,  ,  1  dal  ,  , 

—  ~ — 7U  4"  — W  +  — — y  r — U  W  , 

2c2u  3z  2 c2u  dz 

C0e  ,  1  3cj2  ,  , 

—  ^ T  — — 7  — vw, 

2al  2al  dz 


C0e  ,1  3a^v 

w  +  2ir 


2a2 


w 


— V  +  1 


where  cl,  and  o^,  are  the  Eulerian  velocity  variances  in  x-,y-,  and  z-directions, 
respectively. 


The  meteorological  inputs  of  the  model  are  the  vertical  profiles  of  the  mean  wind  speed  U,  of 
the  Eulerian  velocity  variances  cl,  cl,  and  CJ;,,  and  the  turbulence  kinetic  energy  (TKE)  mean 
dissipation  rate  e.  In  this  report  we  consider  only  a  neutrally  stratified  shear  flow  (i.e.,  Obukhov 
length  L  — >  °°),  and  we  assume  that  the  Eulerian  velocity  variances  in  the  x-  and  y-directions 
(horizontal  velocity  variances)  are  equal  so  cl  =  cl-  The  average  wind  speed  in  the  surface 
layer  above  the  ground  is  described  analytically  by  the  semi-logarithmic  law-of-the  wall: 

(18)  U(z)  =  ^  logffV 

where  w*  is  the  friction  velocity,  k  «  0.4  is  von  Karman’s  constant,  and  zo  is  the  surface 
roughness  length.  Our  parameterizations  for  the  turbulence  velocity  flow  statistics  follow  those 
suggested  by  Rodean  [13]  for  a  neutrally-stratified  boundary-layer  flow: 


(19) 

= 

cl 

(20) 

cl  = 

ul 

“O-sJ'l- 

/  Z\l-5 

(21) 

e  = 

KZ 

(1-0.85-)  , 

where  H  is  the  depth  of  the  (atmospheric)  boundary  layer. 


The  stochastic  differential  equations  (2)  and  (3)  with  the  drift  coefficients  (conditional 
acceleration)  determined  from  Eq.  (17)  are  integrated  forward  in  time  using  a  simple 
first-order  Euler  scheme.  With  this  scheme,  the  update  from  time  t„  to  time  tn+\  (with  time  step 
At  =  tn+ 1  —tn)  for  advancing  the  fluid  particle  position  and  velocity  has  the  form 

(22)  Xn+I  =  X"  +  U"Ar, 

(23)  U^1  -  U"  +  a" At  +  (Coe") 1/2 AW”, 
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where  A W,  is  a  Gaussian  pseudo-random  number  with  mean  (AW,)  —  0  and  covariance 
(AWiAWj)  =  Al 8/ j .  Superscripts  n  and  n  +  1  in  Eq.  (23)  denote  the  values  at  time  tn  and 
tn+\  —tn  +  At,  respectively,  and  a"  is  the  drift  coefficient  vector  evaluated  at  time  tn,  i.e. 
a”  =  a(U”,X”.f„).  Note  that  the  coefficients  in  Eq.  (23)  are  evaluated  at  X”  and  U"  and,  hence, 
the  Euler  approximation  is  explicit.  The  time  step  At  is  chosen  to  be  a  small  fraction  of  the 
Lagrangian  time  scale  Tl  =  2o\./(Cqz)  (more  specifically,  a  time  step  At  —  0.01  T/.  was  found 
to  give  both  an  accurate  and  stable  solution  to  the  stochastic  system  for  fluid  particle 
evolution).  The  value  of  the  universal  constant  Co  possesses  considerable  uncertainty.  The 
existing  literature  seems  to  suggest  that  it  depends  on  the  structure  of  the  turbulent  flow  being 
considered.  As  regards  to  the  choice  of  Co  for  atmospheric  dispersion  applications,  it  appears 
that  a  choice  Co  ~  3.5  provides  good  agreement  with  data  [2,  10]  for  the  one-dimensional  LS 
model  (model  that  neglects  the  shear  stress  of  the  background  flow). 


Simulation  Results _ 

Continuous  point  source 

In  the  first  set  of  numerical  experiments,  we  test  the  code  with  experimental  data  and  compare 
the  estimate  for  the  mean  concentration  obtained  using  the  box  counting  and  the  fast  kernel 
density  estimator  methods  with  measured  data.  To  this  purpose,  concentration  profiles 
measured  in  the  atmospheric  surface  layer  under  near  neutral  conditions  during  the  Project 
Prairie  Grass  experiment  will  be  used  as  a  reference  for  comparison. 

Project  Prairie  Grass  represents  a  benchmark  field  experiment  for  turbulent  dispersion  in  the 
atmospheric  surface  layer.  The  Project  Prairie  Grass  (PPG)  experiment  is  fully  discussed  in  the 
reports  by  Barad  [14]  and  Haugen  [15].  A  continuous  point  source  of  sulfur  dioxide  (SO2)  at  a 
height  of  0.46  m  above  an  extensive  flat  plain  was  used  in  this  experiment.  In  each  run,  the 
time-averaged  (mean)  concentration  was  measured  over  a  sampling  period  of  10  min  using  a 
large  number  of  downwind  detectors  placed  on  circular  arcs  centered  at  the  source.  All  data 
considered  herein  are  from  the  arc  of  detectors  with  a  radius  of  100  m  as  this  was  the  only  arc 
at  which  the  vertical  concentration  profile  was  measured.  Equivalent  two-dimensional 
concentration  profiles  Cy(x,z)  with  x  —  100  m  were  derived  by  performing  a  crosswind 
integration  on  the  concentration;  that  is,  the  crosswind  integrated  concentration  (CWIC)  was 
obtained  from 


(24)  Cy(x,z)  =  J C(x,y,z)dy. 

Sulfur  dioxide  is  almost  a  perfect  conservative  tracer  with  the  uptake  by  the  (grass)  surface 
being  small  (i.e.,  the  surface  is  almost  a  perfect  reflector).  We  therefore  treated  the  surface  as  a 
perfect  reflector  at  the  roughness  length  zo  when  calculating  fluid  element  trajectories.  In  the 
PPG  experiment,  the  roughness  length  of  the  surface  was  estimated  to  be  about  zo  ~  0.006  m. 
In  application  of  the  kernel  estimator  for  the  mean  concentration,  the  concentration  at  a  given 
point  near  the  surface  must  account  explicitly  for  the  reflection  of  particles  from  the  ground 
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surface,  with  the  result  that  the  kernel  estimate  for  C(x,t)  given  by  Eq.  (10)  is  replaced  by 


(25) 


C(xj )  =  Am 


xin)^h)  +K{X~  x(,,)*(0;a) 


where  X^*  =  (xfn\  Y^n\  —  z[n^)  is  the  location  of  a  virtual  fluid  particle  that  is  associated 
with  the  real  fluid  particle  at  the  location  X''7  =  (x[n^ ,  ,  z\n> )  (viz. ,  the  virtual  fluid  particle 

is  the  mirror  image  of  the  real  fluid  particle  with  respect  to  the  ground  surface  at  z  =  0).  More 
specifically,  for  particles  near  the  ground  surface  (at  a  vertical  distance  less  than  the  kernel 
bandwidth  h)  at  X^  =  (x[n\Y[n\z^)  (with  Z'f  <  h),  an  equivalent  virtual  (or,  ghost) 
particle  is  also  introduced  at  X^*  =  (xf”\Yf”\—  z[”^)  to  provide  a  symmetrical  surface 
boundary  condition  at  z  =  0.  It  is  noted  that  for  the  continuous  release  considered  here,  the 
position  of  the  ‘marked’  fluid  particles  are  written  to  the  data  file  every  8T  s  for  the  entire 
simulation  time,  and  this  information  is  used  to  constuct  a  kernel  estimate  for  the  mean 
concentration  exhibited  in  Eq.  (25).  In  particular,  each  ‘marked’  fluid  particle  carries  the 
source  mass  flow  rate  Q/Np  (Q  is  the  mass  flow  rate  of  the  source)  so  the  mass  carried  by  each 
particle  over  the  interval  8T  is  8m  =  Q8T /Np).  For  box -counting,  the  steady-state  mean 
concentration  in  a  small  cell  (Ax  x  Ay  x  Az)  enclosing  the  receptor  point  (x,y,z)  is  estimated  as 
follows: 


(26) 


C(x,y,z) 


Q 

NpAxAyAz 


nP 


where  Q  is  the  source  strength  (kg  s_l)  and  T„  is  the  total  time  the  n-th  particle  remains  in  the 
cell  (residence  time)  and  is  calculated  by  accumulating  the  time  steps  dt  whenever  the  n-th 
particle  resides  in  the  cell. 


Figure  2  shows  the  observed  and  predicted  vertical  profile  of  normalized  cross  wind-integrated 
concentration  utCy / Q  at  a  downwind  fetch  x=  100  m  due  to  a  continuous  source  of  strength 
Q  (g(S02)s-‘)  releasing  tracer  into  a  near-neutral  surface  layer.  With  this  normalization,  the 
normalized  crosswind  integrated  concentration  is  independent  of  w*  (for  a  fixed  zq  and 
experimental  geometry).  The  observations  here  represent  the  average  from  the  nine  PPG  runs 
that  were  the  nearest-neutral  in  atmospheric  stratification  (namely,  the  average  from  PPG  runs 
49,  62,  26,  61,  30,  20,  33,  45,  and  57  with  Obukhov  lengths  L  in  the  range  from  —36  to 
—240  m).  The  predicted  results  for  mean  concentration  were  estimated  from  5,000 
independent  particle  trajectories.  Based  on  only  5,000  independent  trajectories,  the  fast  kernel 
density  estimator  (using  the  quadweight  kernel  function)  for  the  mean  concentration  yielded  a 
good  conformance  with  the  PPG  observations.  The  mean  concentration  estimate  based  on  the 
box-counting  method  (using  bins  Ac  =  2  m  and  Az  =  0.1  m)  also  gave  results  that  closely 
resembled  the  PPG  data,  although  with  only  5,000  independent  trajectories  the  statistical  error 
(or,  variability)  in  the  estimate  is  quite  significant.  In  the  box-counting  method,  the 
root-mean-square  statistical  error  scales  as  Np  . 

Figure  3  is  a  plot  showing  the  mean  concentration  estimate  based  on  box-counting  using 
50,000  independent  particle  trajectories.  The  statistical  error  in  this  box-counting  estimate  for 
the  mean  concentration  is  reduced  considerably  in  comparison  with  that  obtained  using  only 
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5,000  particle  trajectories.  Furthermore,  note  that  the  kernel  estimate  for  the  mean 
concentration  based  on  5,000  independent  particle  trajectories  gives  as  good  a  prediction  for 
the  concentration  as  that  based  on  box  counting  using  50,000  particle  trajectories. 


Instantaneous  point  source 

In  this  subsection,  we  consider  the  tracer  dispersion  resulting  from  an  instantaneous  point 
source  released  in  the  neutrally-stratified,  horizontally  homogeneous,  wall-shear  layer 
(constant  stress  layer).  To  simulate  the  mean  concentration  field  resulting  from  an 
instantaneous  point  source,  thousands  of  particles  are  typically  released  at  the  source  location 
with  their  initial  velocities  equal  to  the  Eulerian  turbulent  velocities  at  this  location,  and  their 
random  trajectories  calculated  using  Eqs.  (2)  and  (3),  with  the  drift  term  a,  (i  =  1,2,3  =  u,  v,  iv) 
specified  by  Eq.  (17).  In  the  following  simulations,  an  instantaneous  point  source  located  at 
x  =  y  =  0  and  z  =  1  m  releases  Q  —  1  kg  of  tracer  into  the  neutrally-stratified  constant  stress 
layer  with  a  roughness  length  zq  —  0.006  m  and  friction  velocity  u*  =  0. 1  ms-1 . 

Figures  4  and  5  show  comparisons  of  the  mean  concentration  estimate  along  the  mean  cloud 
centerline  at  y  =  0  and  at  a  height  z  =  1  m  above  ground  level  obtained  using  the  direct  (naive) 
implementation  of  the  density  kernel  estimator  and  the  fast  kernel  estimator  algorithm 
described  herein.  The  concentration  estimates  were  based  on  50,000  independent  fluid  particle 
trajectories,  and  were  obtained  at  six  different  times  after  the  instantaneous  point  source 
release  (at  t  —  0).  The  density  kernels  used  here  were  the  parabolic  (Epanechnikov)  ( Figure  4) 
and  quadweight  (Figure  5)  kernel  functions  applied  with  a  bandwidth  h  —  2  m.  Note  that  the 
direct  and  fast  algorithms  for  kernel  density  estimation  gave  identical  results  for  the  mean 
concentration  (as  they  should).  However,  for  the  calculation  of  the  mean  concentration  field  at 
513  receptor  locations  for  6  different  travel  times,  the  fast  algorithm  and  direct  method  for 
density  kernel  estimation  required,  respectively,  14  s  and  2800  s  to  complete.  The  calculations 
cited  here  were  carred  out  on  a  PC  with  a  Xeon  processor  (3.2  GHz  clock  speed)  running  Red 
Hat  Linux  9  Operating  System.  For  this  example,  the  fast  algorithm  is  roughly  200  times  faster 
than  the  direct  method  for  density  kernel  estimation,  and  this  speedup  factor  is  expected  to 
increase  for  applications  where  the  number  of  receptor  locations  and/or  independent  fluid 
particle  trajectories  are  greater  than  what  were  used  for  this  simple  example. 

It  was  shown  in  the  previous  subsection  that  the  trajectory-simulation  model  used  here  leads  to 
predictions  of  mean  concentration  that  are  in  very  good  agreement  with  atmospheric 
measurements  of  concentration  obtained  downwind  of  a  continuous  point  source.  There  are  no 
benchmark  concentration  data  for  an  instantaneous  point  source  release  that  can  be  used  to 
compare  with  concentration  estimates  obtained  using  the  fast  kernel  density  estimator.  In 
consequence,  for  the  case  of  an  instantaneous  point  source,  the  “true”  mean  concentration  that 
will  be  used  as  a  reference  for  comparison  will  be  obtained  using  a  very  large  number 
(5  x  106)  particle  trajectories  with  the  “true”  concentration  obtained  using  the  conventional 
box-counting  method  with  small  grid  cells  (cubes)  of  side  length  one  meter. 

Figures  6  and  7  exhibit  streamwise  (x-wise)  cross-sections  of  the  normalized  mean 
concentration  CjQ  estimates  at  6  different  times  after  the  release  obtained  using  the  fast  kernel 
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density  estimator.  These  streamwise  (x-wise)  cross-sections  were  obtained  at 
(y,z)  =  (1  ra.l  m)  and  (y,z)  =  (1  m,3  m),  respectively.  The  estimates  were  obtained  without 
applying  the  virtual  (ghost)  particle  corrections  near  the  ground  surface  as  described  in 
Eq.  (25).  The  estimates  used  the  parabolic  and  quadweight  kernel  functions  with  a  constant 
bandwidth  h  —  2  m.  The  kernel  density  estimates  were  based  on  50,000  independent  particle 
trajectories.  In  spite  of  this,  there  is  generally  a  good  conformance  of  the  quadweight  kernel 
density  estimate  for  the  mean  concentration  with  the  “true”  mean  concentration  based  on 
5  x  106  independent  particle  trajectories  at  both  z  =  1  and  3  m  (cf.  Figures  6  and  7, 
respectively).  The  kernel  density  estimator  produces  relatively  “smooth”  concentration 
estimates  for  even  50,000  particles.  The  parabolic  kernel  density  estimator  underestimated  the 
“true”  mean  concentration  at  z  =  1  m,  although  exhibited  good  conformance  with  the  “true” 
mean  concentration  at  z  =  3  m.  This  is  because  the  shape  of  parabolic  kernel  is  broader  [see 
Figure  1]  than  that  of  the  quadweight  kernel  and  at  z  =  1  m  above  the  ground,  the  parabolic 
kernel  function  with  h  =  2m  has  its  domain  of  support  “truncated”  significantly  by  the  ground 
surface  at  z  =  0  m. 

Figure  8  exhibits  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
C/Q  estimates  at  (y,z)  =  (1  m,  1  m)  and  at  6  different  times  after  the  release,  obtained  using 
the  fast  kernel  density  estimator.  However,  these  results  were  calculated  using  the  kernel 
density  estimator  of  Eq.  (25)  which  accounts  explicitly  for  the  ground  boundary  effects  by 
introducing  for  each  particle  at  a  vertical  distance  less  than  h  (bandwidth)  from  the  ground  an 
associated  virtual  (ghost)  particle  below  the  ground  surface  placed  at  the  mirror  image  location 
relative  to  the  ground  plane  at  z  =  0.  Note  that  the  mean  concentration  estimates  based  on  the 
parabolic  kernel  are  greatly  improved  (cf.  Figure  6)  when  the  ground  boundary  effects  are 
incorporated. 

Figure  9  shows  the  simulation  results  for  streamwise  (x-wise)  cross-sections  of  the  normalized 
concentration  C/Q  estimates  obtained  at  (y,z)  =  (1  m,3  m)  using  a  fast  kernel  density 
estimator  (for  both  parabolic  and  quadweight  kernel  functions  with  h  =  2  m).  These  estimates 
were  obtained  using  only  20,000  independent  particle  trajectories.  In  spite  of  this, 
intercomparison  of  the  kernel  estimates  with  the  “true”  concentration  shows  good  agreement, 
although  not  unexpectedly  the  statistical  variability  in  the  kernel  estimates  are  larger  than  those 
exhibited  in  Figure  7  for  50,000  particles.  Finally,  Figure  10  displays  the  fast  kernel  density 
estimates  of  the  normalized  mean  concentration  C/Q  obtained  in  streamwise  (x-wise) 
cross-sections  through  the  dispersing  cloud  at  (y,z)  —  (1  m,3  m).  This  kernel  density  estimate 
was  obtained  using  only  10,000  independent  particle  trajectories  for  a  parabolic  kernel  with  a 
larger  bandwidth  h  =  6  m.  In  comparison  with  Figure  9,  note  that  increasing  the  bandwith  h 
from  2  m  to  6  m  reduced  the  statistical  error  (variability)  in  the  concentration  estimate,  but 
increased  the  bias.  This  bias  is  largest  near  the  source  release  (cf.  kernel-based  concentration 
estimate  at  t  —  10  s  with  the  “true”  concentration  where  it  is  seen  that  the  peak  mean 
concentration  is  underestimated,  but  the  width  of  the  streamwise  concentration  profile  is 
overestimated).  In  this  example,  a  smoothing  length  h  =  6  m  is  too  large  and  results  in  an 
over-smoothing  of  the  mean  concentration  near  the  source  (viz.,  at  small  travel  times)  where 
changes  in  the  concentration  are  sharp.  However,  further  from  the  source  when  the  mean 
concentration  gradients  are  not  as  large,  the  larger  kernel  bandwidth  leads  to  good  estimates 
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for  the  concentration  with  even  10,000  particle  trajectories. 


Conclusion 


The  kernel  density  estimator  for  the  mean  concentration  field  given  information  embedded  in 
“marked”  fluid  particle  trajectories  provided  an  important  improvement  on  previous  methods 
for  concentration  estimation  based  on  box-counting  techniques  that  were  sensitive  to  the  size 
and  location  of  imaginary  sampling  volumes.  In  spite  of  this,  kernel  estimators  have  rarely 
been  used  in  Lagrangian  particle  modeling  because  the  straightforward  (direct) 
implementation  of  this  method  leads  to  a  computational  effort  that  is  prohibitive  as  the 
computational  work  scales  as  0(N2)  (where  N&Nrfa  Np).  In  this  report,  we  present  and 
implement  an  algorithm  for  the  rapid  evaluation  of  the  kernel  density  estimate  of  the  mean 
concentration  field  that  requires  an  amount  of  computational  effort  proportional  to  N,  and  this 
method  is  insensitive  to  the  precise  distribution  of  the  “marked”  fluid  particles  in  physical 
space  (degree  of  disorder  of  the  fluid  particle  positions  which  result  from  the  stochastic 
trajectory-simulation  model).  In  practice,  speedups  of  two  to  three  orders  of  magnitude  may  be 
expected  in  application  of  the  fast  kernel  density  estimator  compared  to  the  naive  approach, 
depending  on  the  number  of  particles  simulated  and  the  number  of  receptor  locations  at  which 
the  mean  concentration  needs  to  be  determined,  rendering  previously  computationally 
prohibitive  simulations  feasible.  In  this  report,  both  two-  and  three-dimensional  versions  of  the 
fast  kernel  density  estimator  have  been  constructed  and  implemented.  Results  from  a  number 
of  numerical  experiments  have  been  conducted  to  verify  the  algorithm. 

There  are  several  ways  in  which  the  code  for  the  fast  linked  list  based  algorithm  for  the  density 
kernel  estimator  of  the  mean  concentration  field  can  be  either  generalized  and/or  made  more 
efficient.  The  linked  list  data  structure  used  here  is  appropriate  for  the  case  where  the 
bandwidth  of  the  density  kernel  associated  with  each  particle  is  fixed.  In  principle,  the 
bandwidth  h  need  not  be  constant  and  for  many  applications  it  may  be  useful  to  allow  a 
variable  smoothing  length  that  is  dynamically  adapted  so  that  the  number  of  neighboring  fluid 
particles  within  the  domain  of  support  of  a  kernel  function  centered  on  any  particular  particle 
remains  constant.  The  generalization  of  the  proposed  fast  kernel  density  estimator  for  variable 
kernel  bandwidth  h  will  probably  require  that  the  linked  list  data  structure  used  here  be 
replaced  by  a  hierarchy  tree  data  structure  that  can  be  adopted  to  suit  the  needs  of  a  variable 
bandwidth.  In  the  application  of  the  fast  kernel  density  estimator,  the  code  can  be  made  even 
more  efficient  by  pre-computing  the  values  of  the  kernel  function  at  a  large  number  of 
representative  points  and  storing  the  results  in  a  pre-computed  kernel  function  value  look-up 
table.  Finally,  it  may  well  be  worth  exploring  the  parallelization  of  the  code  so  that  it  can  be 
executed  on  computers  with  highly  parallel  architectures  (e.g.,  Beowulf  clusters). 
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Figure  1:  Profile  of  the  non-normalized  Epanechnikov  (parabolic)  (dashed  line)  and  quadweight  (solid  line)  kernels. 
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Figure  2:  Normalized  mean  concentration  estimates  based  on  the  box-counting  and  fast  kernel  density  estimator 
methods  compared  with  observations  obtained  from  the  Project  Praire  Grass  experiment.  The  concentration 
estimates  were  based  on  5,000  independent  particle  trajectories  calculated  using  a  Lagrangian  Stochastic  model. 
The  kernel  function  used  for  the  mean  concentration  estimate  is  the  quadweight  kernel. 
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Figure  3:  Normalized  mean  concentration  estimates  based  on  the  box-counting  and  fast  kernel  density  estimator 
methods  compared  with  observations  obtained  from  the  Project  Praire  Grass  experiment.  The  concentration 
estimates  obtained  using  the  box-counting  and  fast  kernel  density  estimator  were  based,  respectively,  on  50.000 
and  5,000  independent  particle  trajectories  calculated  using  a  Lagrangian  Stochastic  model.  The  kernel  function 
used  for  the  mean  concentration  estimate  is  the  quadweight  kernel. 
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Figure  4:  Simulation  results  for  the  normalized  mean  concentration  along  the  mean  cloud  centerline  aty  =  0  and  at 
height  z  =  1  m  obtained  using  the  direct  and  fast  kernel  density  estimators  for  an  instantaneous  point  source  release 
at  x  =  y  =  0  and  z  =  1  m.  The  mean  concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from 
the  point  source;  namely,  at  (a)  t  =  10,  (b)t  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  —  50,  and  (f)  t  =  60  seconds.  The 
kernel  function  used  for  this  simulation  is  the  parabolic  (Epanechnikov)  kernel  with  a  bandwidth  h  =  2  m.  The  solid 
line  and  open  circles  represent  the  results  obtained  from  the  direct  and  fast  kernel  density  estimators,  respectively. 
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Figure  5:  Simulation  results  for  the  normalized  mean  concentration  C/Q  along  the  mean  cloud  centerline  aty  =  0 
and  at  height  z  =  1  m  obtained  using  the  direct  and  fast  kernel  density  estimators  for  an  instantaneous  point  source 
release  at  x  =  y  =  0  and  z  =  1  m.  The  mean  concentration  is  obtained  at  6  different  times  after  the  tracer  was 
released  from  the  point  source;  namely,  at  (a)  t  =  10,  (b)  t  —  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60 
seconds.  The  kernel  function  used  for  this  simulation  is  the  quadweight  kernel  with  a  bandwidth  h  =  2  m.  The  solid 
line  and  open  circles  represent  the  results  obtained  from  the  direct  and  fast  density  kernel  estimators,  respectively. 
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Figure  6:  Simulation  results  showing  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
C / Q  obtained  at  (y,z)  =  (1  m,  1  m)  for  an  instantaneous  point  source  release  at  x  =  y  =  0  and  z  =  1  m.  The  mean 
concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from  the  point  source;  namely,  at  (a) 
t  =  10,  (b)  t  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60  seconds.  The  normalized  concentration  estimates 
were  obtained  from  50,000  independent  particle  trajectories  using  the  fast  kernel  density  estimator  with  a  parabolic 
(dotted  line)  and  quadweight  (dashed  line)  kernel  functions  with  a  constant  bandwidth  h  =  2  m.  The  solid  line 
corresponds  to  the  "true”  mean  concentration  obtained  from  5  x  106  independent  particle  trajectories  using  a 
box-counting  method.  The  effects  of  the  ground  surface  on  the  kernel  density  estimate  were  not  included  in  this 
simulation. 
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Figure  7:  Simulation  results  showing  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
C / Q  obtained  at  (y,z)  =  (l  m.  3  m)  for  an  instantaneous  point  source  release  at  x  =  y  =  0  and  z  =  1  m.  The  mean 
concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from  the  point  source;  namely,  at  (a) 
t  =  10,  (b)  t  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60  seconds.  The  normalized  concentration  estimates 
were  obtained  from  50,000  independent  particle  trajectories  using  the  fast  kernel  estimator  with  parabolic  (dotted 
line)  and  quadweight  (dashed  line)  kernel  functions  with  a  constant  bandwidth  h  =  2  m.  The  solid  line  corresponds 
to  the  "true"  mean  concentration  obtained  from  5  x  106  independent  particle  trajectories  using  a  box-counting 
method.  The  effects  of  the  ground  surface  on  the  kernel  density  estimate  were  not  included  in  this  simulation. 
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Figure  8:  Simulation  results  showing  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
C/Q  obtained  at  (y,z)  =  (1  m,  1  m)  for  an  instantaneous  point  source  release  atx  =  y=0  and  z  =  1  m.  The  mean 
concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from  the  point  source;  namely,  at  (a) 
t  =  10,  (b)  I  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60  seconds.  The  normalized  concentration  estimates 
were  obtained  from  50,000  independent  particle  trajectories  using  the  fast  kernel  density  estimator  with  parabolic 
(dotted  line)  and  quadweight  (dashed  line)  kernel  functions  with  a  constant  bandwidth  h  —  2  m.  The  solid  line 
corresponds  to  the  “true"  mean  concentration  obtained  from  5  x  106  independent  particle  trajectories  using  a 
box-counting  method.  The  effects  of  the  ground  surface  on  the  kernel  density  estimate  were  included  in  this 
simulation. 
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Figure  9:  Simulation  results  showing  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
CjQ  obtained  at  (y.z)  =  (1  m,3  m)  for  an  instantaneous  point  source  release  atx=y  =  0  and 1  m.  The  mean 
concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from  the  point  source;  namely,  at  (a) 
t  =  10,  (b)  t  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60  seconds.  The  normalized  concentration  estimates 
were  obtained  from  20,000  independent  particle  trajectories  using  the  fast  kernel  density  estimator  with  parabolic 
(dotted  line)  and  quadweight  (dashed  line)  kernel  functions  with  a  constant  bandwidth  h  =  2  m.  The  solid  line 
corresponds  to  the  “true”  mean  concentration  obtained  from  5  x  106  independent  particle  trajectories  using  a 
box-counting  method. 
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Figure  10:  Simulation  results  showing  streamwise  (x-wise)  cross-sections  of  the  normalized  mean  concentration 
C/Q  obtained  at  (y, :)  =  (1  m,3m)  for  an  instantaneous  point  source  release  atx  —  y  =  0  and  z=  1  m.  The  mean 
concentration  is  obtained  at  6  different  times  after  the  tracer  was  released  from  the  point  source ;  namely,  at  (a) 
t  =  10,  (b)  t  =  20,  (c)  t  =  30,  (d)  t  =  40,  (e)  t  =  50,  and  (f)  t  =  60  seconds.  The  normalized  concentration  estimates 
were  obtained  from  10,000  independent  particle  trajectories  using  the  fast  kernel  density  estimator  with  a 
quadweight  (dashed  line)  kernel  function  with  a  constant  bandwidth  h  =  6  m.  The  solid  line  corresponds  to  the  “true” 
mean  concentration  obtained  from  5  x  106  independent  particle  trajectories  using  a  box-counting  method. 
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